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Gait-Based  Human  Recognition  by  Ciassification 
of  Cyciostationary  Processes  on  Noniinear 

Shape  Manifoids 

David  Kaziska  and  Anuj  Srivastava 


We  study  the  problem  of  analyzing  and  classifying  human  gait  by  modeling  it  as  a  stochastic  process  on  a  shape  space.  We  consider  gait 
as  a  evolution  of  human  silhouettes  as  seen  in  video  sequences,  and  focus  on  their  shapes.  More  specifically,  we  define  a  shape  space  of 
planar,  closed  curves  and  model  a  human  gait  as  a  stochastic  process  on  this  space.  Due  to  the  periodic  nature  of  human  walk,  this  process 
is  naturally  constrained  to  be  cyciostationary,  that  is,  its  mean  path  is  assumed  to  be  cyclic.  We  compare  two  subjects  using  a  metric  that 
quantifies  differences  between  average  gait  cycles  of  each  subject.  This  computation  uses  several  tools  from  differential  geometry  of  the 
shape  space,  including  computation  of  geodesics,  estimation  of  means  of  observed  shapes,  interpolation  between  observed  shapes,  and 
temporal  registration  of  two  gait  cycles.  Finally,  we  apply  a  nearest-neighbor  classifier,  using  the  gait  metric,  to  perform  human  recognition, 
and  present  results  from  an  experiment  involving  26  subjects. 

KEY  WORDS:  Biometrics;  Gait  recognition;  Shape  analysis;  Shape  classification;  Statistics  on  shape  manifolds. 


1.  INTRODUCTION 

In  this  work  we  study  the  problem  of  analyzing  videos  of 
humans  with  a  goal  of  recognizing  them  by  analyzing  their 
gait.  In  the  field  of  biometrics,  there  is  a  strong  need  to  recog¬ 
nize  subjects  from  a  distance,  especially  in  noncooperative  en¬ 
vironments.  In  this  situation  images  of  faces,  fingerprints,  or 
irises  may  not  be  available;  a  common  solution  is  to  use  the 
style  of  walking,  called  gait,  to  recognize  people.  Human  gait 
is  valuable  as  a  biometric  because  it  can  be  observed  from  a 
distance  and  requires  no  physical  contact.  Gait  recognition  has 
been  a  problem  of  interest  since  the  seminal  work  of  Niyogi 
and  Adelson  (1994),  where  the  authors  took  advantage  of  the 
periodic  nature  of  gait  and  expressed  an  individual’s  gait  as 
a  combination  of  a  canonical  walk  and  individual  variation. 
In  the  subsequent  literature,  the  approaches  have  been  divided 
into  shape-based  and  motion-based  approaches.  In  the  shape- 
based  approaches,  researchers  use  landmark  representations  or 
binary  images  of  the  subjects  at  various  points  in  the  gait  cycle 
(Liu,  Malave,  and  Sarkar  2004a).  Representations  of  gait  us¬ 
ing  a  mean  shape  on  a  (landmark-based)  shape  manifold  have 
been  used  (Wang,  Ning,  Hu,  and  Tan  2002;  Liu  et  al.  2004b), 
but  without  exploiting  the  gait’s  periodic  nature.  Motion-based 
methods  often  use  parametric  representations  of  certain  bodily 
movements  to  represent  gait.  Some  motion-based  models  com¬ 
bine  parametric  representations  of  arm  and  leg  swing,  stride 
length  and  cadence  and  combine  this  information  with  sta¬ 
tic  data,  such  as  height  (Cuntoor,  Kale,  and  Chellappa  2003; 
BenAbdelkader,  Cutler,  and  Davis  2002,  2004;  Foster,  Nixon, 
and  Prugel-Bennett  2003).  Kale,  Rajagopala,  Cuntoor,  Krueger, 
and  Chellappa  (2004)  constructed  a  hidden  Markov  model  that 
combines  the  shape-  and  motion-based  methods  to  some  extent. 

Our  study  of  gait  recognition  uses  a  shape-based  approach 
based  on  discriminating  shapes  of  silhouettes  of  humans  dur¬ 
ing  walking.  It  is  an  extension  of  a  traditional  statistical  shape 
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analysis,  in  which  single,  individual  shapes  are  compared,  to  a 
comparison  of  sequences  of  shapes.  Shape  analysis  traces  its 
origins  to  the  work  of  D.  G.  Kendall  (1984)  who  used  a  rep¬ 
resentation  of  shape  by  an  indexed  set  of  landmarks  (salient 
points  on  objects)  in  the  plane.  Book-length  treatments  of  this 
approach  include  those  of  Small  (1996)  and  Dry  den  and  Mardia 
(1998).  Landmark-based  representations  have  been  successful 
in  many  applications,  especially  in  physician-assisted  medical 
image  analysis  where  landmarks  are  readily  available,  and  have 
led  to  seminal  advances  in  statistical  analysis  of  shapes  (Dry den 
and  Mardia  1998;  Kent  and  Mardia  2001;  Holboth,  Kent,  and 
Dryden  2002).  Asymptotic  theory  for  this  approach  has  been 
well  developed  by,  among  others,  Bhattacharya  and  Patrange- 
naru  (2002,  2003).  Recently,  it  has  been  extended  to  make  it 
view-independent  by  constructing  a  manifold  of  shapes  that  are 
within  a  projective  transformation  of  one  another  (Mardia  and 
Patrangenaru  2005),  a  thorough  development  of  parametric  dis¬ 
tributions  on  high-dimensional  landmark  shape  spaces  was  pro¬ 
vided  by  Dryden  (2005).  The  study  of  stochastic  processes  on 
shape  manifolds  traces  back  to  work  on  diffusion  on  shape  and 
Brownian  motion  on  shape  manifold  by  D.  G.  Kendall  (1977, 
1984,  1991),  W.  S.  Kendall  (1988,  1990),  and  Le  (1991,  1994). 

Our  approach  to  shape  analysis  essentially  follows  the  frame¬ 
work  of  Klassen,  Srivastava,  Mio,  and  Joshi  (2004)  and  Mio 
and  Srivastava  (2004).  The  basic  idea  is  to  find  a  set  of  all  rel¬ 
evant  closed  curves,  quotient  out  all  shape-preserving  transfor¬ 
mations  from  this  set,  and  use  the  resulting  quotient  space  for 
statistical  shape  analysis.  Those  authors  constructed  shape  man¬ 
ifolds  of  parameterized  curves,  with  individual  parameterized 
curves  represented  by  their  angle  functions  and  speed  functions, 
to  represent  and  analyze  shapes.  Similar  constructions  were 
also  used  by  Younes  (1998, 1999).  Michor  and  Mumford  (2004) 
studied  the  problem  of  comparing  shapes  of  closed  curves  under 
various  metrics.  Using  the  representation  and  metrics  described 
by  Klassen  et  al.  (2004),  Srivastava,  Joshi,  Mio,  and  Liu  (2005) 
presented  techniques  for  clustering,  learning,  and  testing  pla¬ 
nar  shapes.  The  work  of  Klassen  et  al.  (2004)  was  restricted 
to  arc-length  parameterization  of  curves;  Mio  and  Srivastava 
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(2004)  extended  this  by  relaxing  to  permit  variable-speed  para¬ 
meterization,  thus  introducing  a  representation  that  allows  both 
bending  and  stretching  of  shapes  to  match  them  with  one  an¬ 
other.  The  analysis  resulting  from  this  approach  seems  more 
natural  because  interesting  features,  such  as  corners,  are  better 
preserved  while  constructing  statistics  in  this  approach.  Statis¬ 
tical  shape  models  on  this  manifold  have  been  studied  by  Sri- 
vastava,  Jain,  Joshi,  and  Kaziska  (2006). 

Tools  used  in  shape  analysis  often  come  from  differential 
geometry  and  statistics.  To  help  the  reader  with  some  relevant 
concepts  from  differential  geometry,  we  present  a  short  intro¬ 
duction  in  Appendix  A.  We  define  important  quantities,  such 
as  geodesics,  exponential  maps,  group  actions,  and  Karcher 
means,  all  in  a  general  setting.  In  Section  2  we  construct  the 
specific  shape  manifold  of  interest  (i.e.,  the  shape  space  of  elas¬ 
tic  curves),  and  in  Section  3  we  particularize  tools  from  differ¬ 
ential  geometry  to  this  shape  manifold.  In  Section  4  we  present 
our  framework  for  gait  recognition,  in  which  we  explore  the 
concept  of  cyclostationary  processes  on  the  shape  manifold. 
We  define  cyclostationary  processes  and  develop  a  method  for 
their  classification.  We  also  present  important  ingredients  of  our 
framework,  including  automatic  detection  of  gait  cycles,  reg¬ 
istration  of  gait  sequences,  computation  of  a  mean  gait  cycle, 
and  matching  using  geodesic  distances.  In  Section  5  we  illus¬ 
trate  our  experimental  setup  using  a  infrared  (nighttime)  video 
camera  and  present  some  classification  results.  We  successfully 
matched  17  of  26  individuals  for  a  set  of  gait  sequences  that  we 
extracted  from  the  infrared  video. 

2.  SHAPE  SPACE  OF  ELASTIC  CURVES 

Our  approach  to  gait  recognition  is  based  on  modeling  the 
evolution  of  human  silhouettes  as  a  stochastic  process  on  a 
shape  space.  This  space,  comprising  shapes  of  all  simple  closed 
curves  in  is  a  nonlinear  space,  and  tools  from  differen¬ 
tial  geometry  are  needed  to  perform  calculus  on  this  space.  To 
help  the  reader  understand  geometric  terminology  (especially 
a  reader  not  familiar  with  differential  geometry),  we  provide  a 
short  introduction  of  relevant  ideas  from  differential  geometry 
in  Appendix  A. 

The  particular  shape  manifold  that  we  use  in  our  application 
is  the  space  of  elastic  curves  introduced  by  Mio  and  Srivastava 
(2004).  This  space  is  an  extension  of  the  space  constructed  by 
Klassen  et  al.  (2004)  who  used  the  direction  functions  to  rep¬ 
resent  individual  shapes  associated  with  planar,  simple  closed 
curves.  To  develop  this  representation,  start  by  denoting  each 
curve  with  a  differentiable  function  a :  M  ^  which  gives 
parameterized  coordinates.  That  is  to  say  at  a  time  5,  the  vector 
a(s)  e  gives  the  Cartesian  coordinates  of  the  curve.  We  fur¬ 
ther  restrict  our  attention  to  curves  that  are  periodic  with  period 
271  (because  periodicity  makes  curves  closed).  A  closed  curve 
a  is  called  a  simple  closed  curve  if  it  does  not  intersect  itself. 
Our  representation  of  curves  is  in  terms  of  a  pair  of  functions 
(0,  0),  0  :  [0,  Itt]  M,  such  that  at  any  point  5  e  [0,  27r], 

we  have  a'(i')  =  exp(0(5)) exp(y0(5)),  where  j  = 
is  the  instantaneous  speed  of  a  at  5,  and  thus  (j)  is  called  the 
log-speed  function.  The  0{s)  is  the  angle  formed  by  the  vector 
a'(s)  with  the  positive  x-axis,  and,  thus,  0  is  called  the  direction 
function  of  the  curve  a.  For  an  arc-length  parameterized  curve, 
(j)  is  identically  0  but  generally  is  a  nonzero  function. 


Consider  the  space  C  of  all  closed  curves  of  length  In  and 
average  direction  jr  in  given  by 


C 


ds  =  lit, 


1 

—  /  0{s)e‘^^^^  ds  =  it 
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C  is  called  the  preshape  space.  Note  that  the  variability  gener¬ 
ated  by  shape-preserving  transformations  (e.g.,  rotation,  trans¬ 
lation,  and  scale),  are  already  removed  in  this  representation.  In 
other  words,  if  a  curve  a  is  rotated,  translated,  or  scaled,  then 
its  representation  in  C  remains  unchanged.  But,  the  variabil¬ 
ity  resulting  from  different  placements  of  origin  on  the  closed 
curve  a  and  different  reparameterizations  of  [0, 27r]  remain.  In 
other  words,  if  we  place  5  =  0  at  different  points  along  the 
curve,  or  simply  reparameterize  the  interval  [0, 27r],  then  we 
will  get  a  new  pair  (<^,  0)  e  C,  but  the  shape  of  the  underlying 
a  will  remain  unchanged.  We  identify  all  elements  of  C  that 
have  the  same  shape  using  algebraic  equivalence  as  follows. 
The  variability  resulting  from  different  placements  of  origin  can 
be  modeled  as  a  group  action  of  the  unit  circle,  on  C,  and 
that  reparameterization  results  from  the  group  action  of  V,  the 
set  of  all  automorphisms  {y  :  [0,  27r]  [0, 27r]},  on  C.  For  a 

definition  of  a  group  action  on  manifolds,  refer  to  Appendix  A. 
The  action  of  S*  on  C  is  given  by  the  following:  For  a  5o  e 
the  curve  (f(s),0(s))  becomes 


50  •  (fis),  0(s))  =  ((p(s-  5o),  0(s  -  5o)). 


An  automorphism  y  changes  the  representation  of  a  curve  ac¬ 
cording  to 


{f,0)oY  =  {foY  -l-logy',  0  o  y). 


All  of  the  curves  generated  by  changing  the  origin  or  repara¬ 
meterizing  the  same  curve  are  considered  to  be  of  the  same 
shape.  Therefore,  the  shape  space  is  defined  to  be  a  quotient 
space  S  =  C/(S*  x  V).  In  S,  each  point  represents  a  distinct 
shape,  and  S  becomes  the  space  in  which  statistical  analysis  of 
shapes  is  to  be  performed.  In  the  next  section  we  particularize 
tools  from  differential  geometry  to  this  shape  space  S  for  use 
in  our  application  on  gait  recognition. 


3.  TOOLS  FOR  ANALYSIS  ON  SHAPE  SPACE  5 

Our  representation  of  a  human  gait  consists  of  a  temporal  se¬ 
quence  of  simple,  closed  curves,  which  are  the  silhouettes  of 
individuals  in  video  sequences.  Focusing  on  the  shapes  of  these 
silhouettes,  we  consider  them  to  be  points  on  the  shape  space  S. 
Our  method  of  gait  recognition  exploits  the  differential  geom¬ 
etry  of  S  various  ways,  including  interpolating  gait  sequences, 
registering  gait  sequences,  computing  a  mean  gait  cycle,  and 
performing  nearest-neighbor  matching  for  gait  recognition.  In 
all  of  these  tasks,  our  understanding  of  the  differential  geom¬ 
etry  of  S  plays  a  central  role.  We  note  that  the  our  preshape 
space  C  is  a  complete  Riemannian  manifold  and,  consequently, 
between  any  two  points  of  this  space,  a  geodesic  can  be  de¬ 
fined  and  computed.  Furthermore,  because  the  shape  space  S  is 
a  quotient  space  of  C,  geodesics  in  S  are  computed  as  shortest 
geodesics  connecting  equivalence  classes  in  C. 
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Figure  1.  Ten  iterations  of  the  shooting  method  for  computing  geodesics.  In  each  iteration,  the  shooting  direction  is  corrected  to  get  closer  to 
the  target  point. 


3.1  Geodesic  Paths  on  Shape  Space 

An  efficient  technique  for  quantifying  shape  differences  is  to 
compute  geodesic  paths  in  S  connecting  shapes,  and  then  use 
their  lengths  to  quantify  shape  differences.  Because  5  is  a  quo¬ 
tient  space  of  C,  the  geodesic  length  between  any  two  shapes  in 
S  is  given  by 

dsii(t>i,0i),{(l)2,92))=  min  dc(((^i,  0i),  50  •  (02,  ^2)  o  y). 

Here  dc  and  ds  are  the  distances  between  points  on  C  and  S. 
This  equation  states  that  we  fix  one  shape  (0i ,  9\ )  and  seek  the 
reparameterization  of  the  other  shape  (02,  92)  that  minimizes 
the  geodesic  distance  from  (0i ,  0i)  in  C.  The  minimization  over 
50  can  be  easily  performed  using  an  exhaustive  search  over  S^, 
but  the  search  for  optimal  y  deserves  a  closer  look.  We  seek  a 
reparameterization  function  y  sV  of  (02,  ^2)  such  that  it  min¬ 
imizes  the  matching  cost 

H(y)  =  r  (A  II  (01  (5),  01  (5))  -  (02(5),  02(5))  oyf 
Jo 

+  a-X)\y'(s)\^)ds, 


0  <  y  <  I,  and  for  a  fixed  A  >  0.  The  minimization  is  per¬ 
formed  using  the  dynamic  programming  algorithm,  which  is 
briefly  outlined  in  Appendix  B.  Figure  2  shows  some  examples 
of  this  matching.  Note  that  in  the  figure,  the  matching  process 
works  well  whether  the  legs  are  apart  or  together,  hands  are 
visible  or  not  visible,  and  so  on.  This  estimation  of  y  is  also 
called  the  registration  of  one  curve  to  another.  This  registra¬ 
tion  of  curves  is  similar  in  concept  to  the  registration  of  gait 
sequences  that  we  introduce  in  the  next  section,  although  the 
two  cases  differ  in  the  actual  quantities  being  registered. 

Once  the  two  curves  are  registered,  they  can  be  treated  as 
elements  of  C,  and  the  computation  of  a  geodesic  is  based  on 
a  shooting  method  described  in  Appendix  A.  To  illustrate  the 
resulting  geodesics,  Figure  3  shows  three  examples  of  geo¬ 
desic  paths  between  shapes  of  human  silhouettes.  In  each  row, 
the  shapes  denote  equally  spaced  points  along  a  geodesic  con¬ 
necting  the  starting  and  ending  shapes.  For  any  two  shapes 
(01 , 0i)  and  (02,  02),  the  length  of  geodesic  path  between  them 
forms  a  natural  tool  for  comparing  them  and  is  denoted  by 
dsii4>i,  0i),  (02,  02))-  Also,  we  use  the  function  'h(t)  to  denote 


Figure  2.  Shape  matching.  Each  pair  shows  a  y  that  minimizes  the  matching  cost  between  the  two  shapes;  the  lines  show  the  corresponding 
points  between  the  two  curves. 
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Figure  3.  Examples  of  geodesic  paths  between  human  shapes  in  S.  In  each  row,  the  first  and  the  last  shapes  are  given,  and  the  intermediate 
shapes  depict  points  along  geodesics  connecting  these  given  shapes. 


the  geodesic  path,  so  that  ^(0)  =  ((^i,  9i)  and  =  (02>  Oi)- 
Here  'h(O),  the  initial  velocity  vector  along  this  path,  is  useful 
in  computing  mean  shapes. 

3.2  Piecewise-Geodesic  Interpolations  Between  Shapes 

Our  framework  considers  the  underlying  gait  sequence  as  a 
stochastic  process  on  S.  In  practice,  one  has  only  a  finite,  dis¬ 
crete  set  of  observations  along  a  sample  path  when  the  human 
walk  is  captured  using  a  video  camera.  For  the  purpose  of  align¬ 
ing  (or  registering)  different  sequences,  as  described  later,  it 
will  be  useful  to  fill  in  between  the  observed  shapes  to  result 
in  dense  time  samples.  We  use  a  piecewise-geodesic  interpo¬ 
lation  to  fill  in  between  the  observed  shapes.  In  other  words, 
we  connect  each  successive  pair  of  observed  shapes  with  geo¬ 
desic  paths  on  S,  and  in  this  way  we  can  obtain  an  arbitrarily 
dense  time  sampling.  Interpolation  is  an  important  tool  in  our 
framework  because  it  allows  us  to  analyze  a  gait  process  at  arbi¬ 
trary  time  resolutions.  To  compare  two  gait  sequences,  we  use 
a  metric,  given  in  (3),  that  requires  samples  of  sequences  at  si¬ 
multaneous,  discrete  points  in  time.  When  the  time  points  do 
not  correspond  to  the  times  at  which  shapes  are  observed,  we 
will  need  to  interpolate  on  S  to  estimate  the  silhouette  at  those 
times.  To  guard  against  bias  that  can  be  introduced  in  this  inter¬ 
polation,  we  performed  a  limited  amount  of  interpolation.  We 


performed  only  the  amount  of  interpolation  needed  to  standard¬ 
ize  the  number  of  silhouettes  in  a  gait  cycle. 

Suppose  that  we  wish  to  interpolate  between  two  con¬ 
secutive  observed  silhouettes  having  representations  (</>/,  0,) 
and  (</>i+i,  ft+i).  Using  the  shooting  method  for  computing 
geodesics,  we  compute  a  geodesic  between  the  two  points. 
The  resulting  geodesic  ^  from  will  reach  the  target 

shape  (</>/+!,  0i+i)  in  unit  time,  'k(O)  =  {(pi,9i)  and  vi/(i)  = 
We  can  then  evaluate  'k(t)  for  any  t  €  (0,  1)  to 
continuously  interpolate  between  the  two  observed  shapes. 

To  demonstrate  this  idea  experimentally,  we  present  an  ex¬ 
ample  in  Figure  4.  We  start  with  a  gait  sequence  consisting  of 
six  silhouettes,  which  represents  a  half-cycle  of  gait  for  this 
individual.  To  illustrate  geodesic  interpolation,  we  first  drop 
the  third  shape  in  this  sequence  and  use  the  second  and  fourth 
shapes  to  estimate  it  using  interpolation.  The  first  contour  in 
the  second  row  is  this  interpolated  shape.  Note  its  closeness 
to  the  actual  observed  third  shape  in  the  original  sequence.  The 
same  experiment  is  repeated  for  the  fourth  shape;  it  is  estimated 
by  interpolating  between  the  third  and  fifth  shapes.  The  esti¬ 
mated  shape  in  shown  as  the  second  shape  in  the  second  row. 
Because  alternate  techniques  in  the  literature  involve  binary  im¬ 
ages,  we  performed  an  interpolation  on  the  binary  images  pro¬ 
duced  from  the  second  and  fourth  silhouettes;  the  resulting  in- 


Sequeiice  of  Observed  Shapes 


Interpolated  Sequence 


Figure  4.  Illustration  of  gait  cycles.  The  sequence  from  legs  together  to  right  leg  forward  to  legs  together  is  the  first  half-cycle.  The  sequence 
from  legs  together  to  left  leg  forward  to  legs  together  is  the  second  half-cycle. 
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Sequence  of  Observed  Shapes 


Interpolated  Sequence 


Figure  5.  Detection  of  gait  cycles.  On  the  left  is  an  ordered  se¬ 
quence  of  silhouettes.  On  the  right  is  a  plot  of  the  geodesic  distance 
from  the  first  shape  to  each  subsequent  shape.  The  peaks  at  shapes  5, 
12,  18,  and  25  correspond  to  the  beginnings  and  ends  of  cycles. 

terpolated  gray-scale  image,  is  shown  in  the  third  row  of  Fig¬ 
ure  4.  Similarly,  we  interpolated  between  the  third  and  fifth  sil¬ 
houettes  and  compared  both  methods  with  the  observed  fourth 
silhouette. 

3.3  Computation  of  Mean  Shapes 

In  the  proposed  statistical  framework  for  gait  representation 
and  analysis,  another  important  ingredient  is  computation  of 
the  statistical  mean  of  observed  shapes.  In  other  words,  given 
a  finite  collection  of  points  on  S,  we  would  like  to  define  and 
compute  a  quantity  that  represents  the  central  tendency  of  that 
dataset.  For  this  purpose,  we  use  the  notion  of  a  Karcher  mean, 
which  is  essentially  the  centroid  of  the  sampled  data  points  on 
S  using  the  geodesic  distance  d^.  This  quantity  has  also  been 
called  the  Frechet  mean  or  the  intrinsic  mean  in  the  litera¬ 
ture.  Given  a  set  of  sample  shapes  ((pi,  9i), . . . ,  (cpn,  On)  €  S, 
the  sample  Karcher  variance  is  a  function  of  ((p,0)  €  S  and 
is  given  by  V((p,  0)  =  0),  {(pt,  Ok))^-  The  Karcher 

mean  set  of  {(pi,0i), . . . ,  {cpn,  9„)  e  S  is  the  set  of  minimiz- 
ers  of  V{9).  Several  authors,  including  Klassen  et  al.  (2004), 
have  used  an  iterative  algorithm  for  computation  of  the  sample 
Karcher  mean. 

4.  HUMAN  IDENTIFICATION  USING  GAIT  ANALYSIS 

In  this  section  we  present  a  statistical  framework  for  per¬ 
forming  gait  analysis  and  its  use  in  human  recognition.  Human 
gait  is  modeled  as  a  stochastic  process  on  the  shape  space  S. 
With  this  formulation,  we  are  representing  silhouettes  of  peo¬ 
ple  walking  using  simple  closed  curves.  We  recognize  that  there 
can  be  times  when  armswing  or  leg  movements  prevent  a  sil¬ 
houette  from  being  a  simple  closed  curve.  We  did  not  observe 
such  cases  in  practice,  but  in  this  framework  we  would  repre¬ 
sent  this  situation  by  the  simple  closed  curve  obtained  from  the 
outermost  contour. 

Because  human  walk  is  quite  naturally  a  periodic  process,  we 
restrict  ourselves  to  a  cyclostationary  process,  defined  next. 

Definition  1  (Cylostationary).  A  stochastic  process  is  called 
cyclostationary  with  a  period  r  if  the  joint  probability  distribu¬ 
tion  of  the  random  variables  X{ti),  X  (t2), . . X  fin)  is  same  as 
that  of  Xfi\  +  x),Xfi2  +  T), ...  ,Xfin  +  t),  for  all  ti ,  t2,  ■  •  ■ ,  G 
and  for  all  n.  In  particular,  the  random  quantities  Xfi)  and 
Xfi  +  t)  have  the  same  probability  distribution. 
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Sequence  2  -  Before  Registration 


Sequence  2  -  After  Registration 


Figure  6.  Observed  gait  cycles  for  two  different  people.  The  cor¬ 
respondence  problem  determines  which  shapes  in  the  first  and  second 
sequences  are  to  be  compared. 

Previous  applications  of  cyclostationary  processes  have  in¬ 
cluded  the  fields  of  signal  processing  and  meteorological  sci¬ 
ence.  These  applications  have  mostly  involved  low-dimensional 
signals  in  Euclidean  spaces.  To  the  best  of  our  knowledge, 
ours  is  the  first  use  of  the  cyclostationary  structure  on  infinite¬ 
dimensional,  nonlinear  manifolds. 

We  model  the  shape  process  generated  by  temporal  evolution 
of  human  silhouettes  as  a  cyclostationary  process  on  S.  For  any 
cyclostationary  process,  the  notion  of  a  cycle  is  central  to  its 
analysis.  For  a  gait  process,  we  consider  afull  cycle  as  the  pe¬ 
riod  starting  from  when  legs  and  hands  are  all  together  to  the 
time  of  return  to  a  similar  state,  as  shown  in  Figure  7.  The  top 
row  of  Figure  7  shows  the  first  half-cycle  where  the  left  foot 
goes  forward  and  the  right  foot  catches  up,  and  the  bottom  row 
shows  the  second  half-cycle  where  the  right  foot  moves  first.  In 
our  illustrations  in  this  section  and  for  our  classification  results 
in  the  next  section,  we  use  half-cycles.  We  assume  that  gait  se¬ 
quence  associated  with  a  person  is  a  cyclostationary  process  on 
S.  The  duration  of  a  cycle  corresponds  to  the  period  r  of  the 
process.  Like  any  cyclostationary  process,  it  suffices  to  study  a 
gait  sequence  within  the  period  [0,  t].  As  a  result  of  the  cyclo¬ 
stationary  nature  of  gait,  in  our  notation  we  generally  consider 
the  time  t  to  be  modulo  r,  so  that  t  e  [0,  r).  Given  two  stochas¬ 
tic  processes,  our  main  goal  is  to  quantify  differences  between 
them.  Let  X  fi)  and  Y  fi)  be  two  gait  processes  on  the  shape 
space  S,  with  periods  and  Xy.  We  seek  a  metric  dpfiX,  Y) 
with  certain  desired  properties.  There  are  two  possibilities  for 
such  a  metric:  mean  of  squared  distances  and  squared  difference 
between  average  paths. 

Mean  of  Squared  Distances.  The  first  idea  is  to  compute 
the  expected  value  of  the  squared  distance  between  the  sample 
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(13)  (32)  (51) 


(52)  (67)  (90) 


Figure  7.  Demonstration  of  interpolation  between  points  in  S.  The  first  row  shows  the  original  gait  sequence  of  six  shapes.  The  second  row 
shows  results  of  interpolation  on  the  shape  space  for  the  third  and  fourth  shapes.  The  interpolated  shapes  are  drawn  directly  below  the  original 
third  and  fourth  shapes  for  comparison.  In  the  third  row  we  interpolate  the  corresponding  binary  images,  producing  the  gray-scale  images  that 
correspond  to  the  same  shapes. 


paths  of  X  and  Y, 
dpiX,  Y) 


=  min  E 


(^l\s{X{t),Y{K  +  git))fdt 


1/2 


(1) 


where  E  denotes  the  expected  value,  ds(-,  ■)  denotes  the  geo¬ 
desic  length  metric  defined  on  S,  g  denotes  a  smooth  mapping 
between  [0,  t^]  and  [0,  ty],  and  k  denotes  a  possible  relative 
time  shift  between  the  two  observed  gait  sequences.  The  map¬ 
ping  g  is  needed  to  register  shapes  along  the  two  sequences 
before  comparison. 


Squared  Distance  Between  Average  Paths.  The  second  idea 
is  to  use  squared  distance  between  the  mean  cycles, 

dp(X,Y)=  min^  ds{E[X{t)], 

9  \l/2 

E[YiK  +  g(t))]ydt]  ,  (2) 

where  expectations  ^[^(t)]  and  £[T(t)]  are  computed  on  the 
shape  space  S. 

We  note  that  these  functionals  may  not  be  symmetric;  how¬ 
ever,  in  either  case  a  true  metric  can  be  constructed  using 
d(X,  Y)  +  d(Y,  X).  Although  both  distances  are  relevant,  the 
latter  is  more  efficient  from  a  computational  standpoint.  The 
first  definition  requires  computing  distances  between  several 
random  realizations  of  gait  cycles  of  X  and  Y.  In  other  words, 
we  need  to  compute  distances  between  all  training  and  test  cy¬ 
cles  associated  with  X  and  Y.  In  contrast,  the  second  definition 
requires  computing  an  average  gait  cycle  each  for  X  and  Y, 
then  computing  a  distance  between  the  average  cycles.  In  the 
context  of  human  recognition  using  a  distance-based  classifier, 
X  denotes  the  gait  process  of  a  known  person  (training  data) 
and  Y  denotes  that  of  an  unknown  person  (test  data).  Therefore, 
once  the  average  gait  cycle  for  the  known  person  is  computed 
during  training,  it  can  be  simply  stored  for  future  classification 
purposes.  It  can  be  viewed  as  a  template  gait  cycle  for  this  per¬ 
son;  this  is  similar  to  the  idea  of  deformable  template  theory  for 
classification. 


Consequently,  we  choose  the  second  metric  and  use  its  dis¬ 
crete  form. 


dp{X,Y)=  min  ['^d{X{t),Y{K  +  gf))]" 


1/2 


(3) 


u=l 


where  X  and  Y  are  sample  Karcher  means  of  the  processes  X 
and  Y .  With  a  slight  abuse  of  notation,  we  continue  to  use  Xx 
and  Xy  as  lengths  of  gait  cycles  in  the  discrete  case  as  well. 
Now  g  :  {1, 2, . . . ,  T^}  1-^  {1, 2, . . . ,  Ty}  is  a  mapping  that  regis¬ 
ters  points  across  the  two  sequences.  Estimation  of  k  and  g  is 
discussed  in  the  next  two  sections. 


4.1  Automatic  Detection  of  Cycles 

The  first  issue  in  computing  the  distance  given  in  (3)  is  the 
extraction  of  gait  cycles  for  each  person.  We  typically  have 
large  sequences  of  human  silhouettes  generated  from  videos, 
and  we  need  to  extract  a  few  gait  cycles  for  analysis.  Because 
the  process  is  repetitive,  the  geometry  of  the  shape  space  pro¬ 
vides  a  method  for  automatically  detecting  cycles,  because  sim¬ 
ilar  shapes  should  occur  at  the  same  point  in  the  cycle.  In  our 
application  to  human  gait,  we  noticed  that  the  silhouettes  with 
the  arms  and  legs  together  are  far  away,  in  terms  of  geodesic  dis¬ 
tance,  from  the  points  at  which  the  limbs  are  extended.  To  iden¬ 
tify  the  beginnings  and  ends  of  cycles  in  a  sequence  of  observed 
shapes,  we  begin  with  a  silhouette  with  the  limbs  extended,  then 
compute  the  geodesic  distances  from  that  first  shape  to  all  of 
the  following  shapes.  Because  the  shapes  with  arms  and  legs 
together  are  far  from  shapes  with  limbs  extended,  then  the  dis¬ 
tance  that  we  compute  shows  peaks  at  the  shapes  with  the  limbs 
together.  We  detect  the  beginnings  and  ends  of  cycles  by  look¬ 
ing  for  these  peaks.  An  example  of  this  is  shown  in  Figure  8. 

This  automated  detection  of  cycles  obviates  the  need  for  min¬ 
imizing  over  K  in  (3).  The  resulting  metric  is  now 


dliX,  Y)  =  mml^d{X(t),  Y(g{t))y 


1/2 


(4) 


u=l 


JASA  jasa  v. 2007/01/31  Prn:30/05/2007;  15:52 


F: jasaap06263rl . tex;  (Aust)  p.  7 


Kaziska  and  Srivastava:  Gait-Based  Human  Recognition  by  Ciassification 

il  il  il  iJ  H  s 

H  k  A  A  I  (! 

il  il  i  il  i 

(!  H  /I  il  i  i 

S  &  k  A  A  ii 


7 


Figure  8.  Example  of  interpolation  by  linear  time  scaling.  The  first  row  shows  a  sequence  of  eight  silhouettes  taken  from  a  video  clip.  The 
second  row  shows  the  interpolation  and  the  uniform  resampling  of  this  sequence  to  result  in  a  sequence  of  length  six. 


4.2  Registration  of  Gait  Cycles 

Once  we  have  a  method  for  automatically  extracting  cycles, 
the  next  important  problem  in  comparison  and  recognition  of 
processes  is  the  registration  of  points  along  two  processes.  To 
formalize,  let  Xx  and  Xy  be  the  periods  of  two  processes,  X  and 
Y,  and  let  g  :  [0,  r^]  [0,  Ty]  be  a  map  that  is  invertible,  with 

both  g  and  g~^  having  piecewise  continuous  derivatives.  Our 
goal  is  to  find  g’s  that  minimize  energy  functions  such  as  those 
given  in  (l)-(4).  In  a  discrete  implementation,  the  registration 
of  any  two  discrete  cycles  amounts  to  solving  the  equation 

g  =  argmm'^d{X(t),Y(g(t))f. 
s  /=i 

To  explain  this  problem  further  in  the  context  of  human  gait 
cycles,  given  samples  of  shapes  along  two  observed  walks,  we 
need  to  determine  which  shapes  to  compare.  Even  though  the 
shapes  form  an  ordered  sequence,  there  may  be  a  time  scal¬ 
ing,  time  warping,  and/or  time  shifting  between  the  two  se¬ 
quences.  Figure  9  shows  examples  of  observed  gait  cycles  for 
two  people.  One  cycle  contains  10  shapes,  and  the  other  con¬ 
tains  8  shapes.  To  compare  these  two  sequences,  we  need  to 
determine  which  shape  in  the  second  sequence  corresponds  to 
a  given  shape  in  the  first  sequence.  The  shape  that  we  need  may 
be  an  observed  shape  from  the  second  sequence  or  may  be  a 
shape  occurring  between  two  shapes  from  the  first  sequence. 
In  the  latter  case,  we  may  need  to  interpolate  to  estimate  the 
shape  that  we  need.  Thus  two  issues  necessitate  the  registration 
of  processes: 


Figure  9.  Example  of  registration  by  DTW.  The  first  row  shows  a 
sequence  of  10  silhouettes  from  a  training  sequence.  The  second  row, 
shows  a  test  sequence  that  we  wish  to  register.  The  third  ???  row  shows 
the  sequence  after  registration  and  the  function  (p,  found  by  dynamic 
programming,  used  in  the  registration. 


•  Registration  of  two  mean  cycles  to  compare  them.  In  com¬ 
puting  the  metric  in  (4),  we  compare  two  mean  cycles  at 
points  in  time.  It  may  occur  that  process  speeds  within  the 
cycles  vary  between  the  mean  sequences  that  we  compare. 
If  so,  then  registering  the  mean  cycles  may  provide  a  more 
accurate  match. 

•  Registration  of  processes  to  compute  a  mean  cycle.  Reg¬ 
istration  of  observed  processes  is  necessary  to  compute  a 
mean  cycle  for  a  particular  observed  process.  The  cycles 
that  we  observe  in  an  observed  process  may  have  differ¬ 
ent  lengths  if  the  process  changes  speed.  If  the  process 
changes  speed  within  a  cycle,  then  we  may  wish  to  cor¬ 
rect  for  the  speed  change.  For  both  of  these  purposes,  we 
must  register  the  observed  sequences  within  the  sequence 
before  estimating  the  mean  cycle. 

We  achieve  registration  using  one  of  two  techniques:  linear  time 
scaling  or  time  warping  using  dynamic  programming. 

Linear  Time  Scaling.  The  simplest  idea  is  to  consider  g 
simply  as  a  linear  time  scaling,  g(t)  =  pt,  where  /I  >  0  is  a 
scalar.  When  the  endpoints  of  the  two  cycles  are  known,  p  is 
simply  the  ratio  XyjXx-  The  underlying  assumption  here  is  that 
speeds  of  two  processes  X  and  Y  are  constant  during  their  ob¬ 
served  cycles,  and  we  can  match  them  by  simply  rescaling  time- 
wise.  We  illustrate  this  idea  with  an  example  in  Figure  5.  Con¬ 
sider  that  a  half-cycle,  Y ,  shown  in  the  top  row,  has  come  from 
an  observed  sequence,  and  we  wish  to  register  this  sequence  in 
time  to  another  half-cycle  of  length  6  X  (not  shown).  We  ob¬ 
tained  the  silhouettes  in  the  top  row  at  eight  uniformly  spaced 
points  in  time,  which  we  denote  by  t  =  0-7.  Using  p  =  5/1, 
we  can  register  the  first  half-cycle  to  the  second  half-cycle, 
but  we  need  silhouettes  at  the  six  uniformly  spaced  points, 
t  =  0,  -y,  y,  y,7,  for  process  Y.  To  estimate  the  shapes  oc¬ 

curring  at  noninteger  units  of  time,  we  use  piecewise  geodesic 
interpolation  on  Y ;  for  example,  to  get  a  point  at  time  1  =  5, 
we  travel  for  time  .4  along  the  geodesic  path  from  the  shape  at 
r  =  1  to  the  shape  at  t  =  2.  The  remaining  shapes  in  the  second 
row  of  Figure  5  are  computed  similarly. 

Time  Warping  Using  Dynamic  Programming.  In  cases 
where  speeds  of  processes  vary  within  their  cycles,  the  dynamic 
time  warping  (DTW)  approach  (Bellman  2003)  provides  an  al¬ 
ternative  to  linear  interpolation  that  may  provide  a  more  accu¬ 
rate  match.  Given  processes  X  and  Y  on  intervals  [0,  xx]  and 
[0,  Xy],  seeks  a  diffeomorphism,  g  :  [0,  Xx]  [0,  xy],  such  that 
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Interpolated  Cycle  One  Interpolated  Cycde  Two 


Interpolated  Cycle  Three  Mean  Cycle 


Figure  10.  Computation  of  a  mean  cycle.  The  first  three  consist  of  gait  cycles  registered  using  linear  interpolation.  The  mean  cycle  is  in  the 
fourth  set.  Each  shape  in  the  fourth  row  is  the  Karcher  mean  of  the  three  corresponding  shapes. 


the  metric  in  (5)  is  minimized.  Rather  than  a  linear  form  of  the 
function  g,  we  generalize  to  a  function  that  can  have  the  effect 
of  speeding  up  or  slowing  down  an  observed  cycle,  as  neces¬ 
sary.  We  may  need  to  do  this  before  computing  a  mean  cycle 
if  the  subject  changes  speeds  during  a  cycle.  This  technique 
also  may  improve  our  ability  two  match  two  observed  mean 
cycles.  Using  techniques  described  earlier,  it  is  often  possible 
to  directly  detect  the  start  and  end  of  a  shape  gait  cycle  in  a 
given  sequence.  Therefore,  once  the  cycles  have  been  detected, 
the  need  to  estimate  k  is  obviated,  and  we  can  directly  estimate 
g  using  dynamic  programming.  Several  authors  have  used  this 
DTW  idea  in  different  contexts.  This  solves  for 

g  =  argmin'^d{X(t),Y(g(t))fdt.  (5) 

^  ^=1 

Figure  6  illustrates  DTW  showing  one  half-cycle  from  two  dif¬ 
ferent  gait  sequences  X  and  Y  for  which  we  wish  to  register  the 
second  sequence  to  the  first.  We  show  the  registered  sequence 
and  note  that  it  provides  a  visually  better  match  in  time  for  the 
first  sequence,  and  we  plot  the  function  g  that  we  used  to  per¬ 
form  the  registration. 

4.3  Computation  of  Mean  Gait  Cycles 

To  use  (4)  for  comparing  and  classification  of  processes,  we 
need  to  estimate  the  Karcher  mean  shape,  E[X{t)\,  for  rele¬ 
vant  times  in  a  cycle.  Assuming  that  we  have  observed  multiple 
observations  cycle  of  each  process,  the  first  task  is  to  register 
the  shapes  across  cycles,  as  described  earlier,  and  then  compute 
the  means  of  the  corresponding  shapes  in  S.  The  mean  shape 
E[X (?)]  is  defined  as  the  sample  Karcher  mean  shape  at  time  t 
and  computed  as  explained  in  Section  3.3.  Figure  10  shows  an 
example  of  a  calculated  mean  gait  cycle. 

5.  EXPERIMENTAL  RESULTS 
5.1  Experimental  Setup 

As  an  application  of  gait  analysis,  we  used  a  night-vision  or 
infrared,  (IR)  video  camera  to  observe  human  gait.  In  addition 


to  collecting  data  in  noncooperating  environments,  and  from  a 
distance,  an  IR  camera  provides  the  added  benefit  of  working 
in  dark.  This  is  especially  well  suited  to  surveillance  and  ur¬ 
ban  battlefield  scenarios,  where  a  major  portion  of  imaging  is 
done  using  IR  sensors.  Figure  1 1  shows  an  example  of  video 
sequence  obtained  using  an  IR  camera.  The  task  of  extracting 
human  silhouettes  from  video  sequences  was  semiautomated; 
that  is,  it  required  some  human  intervention  in  addition  to  auto¬ 
mated  programs.  We  do  not  describe  this  process  of  silhouette 
extraction  and  assume  that  the  data  are  available  as  sequences 
of  silhouettes  for  each  human  subject.  Our  experimental  results 
are  based  on  a  collection  of  IR  video  clips  of  26  student  and 
faculty  volunteers  from  Florida  State  University.  We  collected 
at  least  two  clips  of  each  person  and  formed  disjoint  training 
and  tests.  We  performed  a  gait-matching  experiment  following 
these  steps: 

•  For  each  of  training  and  test  sequence,  we  extracted  three 
half-cycles,  performed  registration  using  linear  time  scal¬ 
ing,  and  then  computed  an  average  gait  cycle. 

•  For  each  test  sequence,  we  computed  the  metric  in  (3)  for 
each  training  sequence  and  sought  the  nearest  match. 

We  observed  anecdotally  that  the  silhouettes  in  individuals’ 
gaits  were  different  depending  on  which  leg  was  leading,  usu¬ 
ally  due  to  different  armswings.  As  a  result,  the  half-cycles  that 
we  used  were  from  a  single  side  for  each  individual. 

5.2  Classification  Results 

The  classification  results  are  summarized  in  Table  1.  Un¬ 
der  the  nearest-neighbor  (NN)  criterion,  we  obtain  a  successful 
match  for  17  of  the  26  test  sequences.  For  three-NN  classifiers, 
where  the  correct  class  in  the  training  set  must  be  among  the 
three  neighbors  of  the  test  cycle,  the  success  rate  was  21  out  of 
26.  An  example  of  a  correct  match  is  shown  in  the  top  row  of 
Figure  12,  whereas  an  incorrect  match  is  shown  in  the  bottom 
row. 

We  implemented  a  few  current  approaches  to  compare  them 
with  our  results.  First,  we  implemented  a  simple  method,  called 


Figure  11.  A  small  portion  of  a  sample  IR  video  sequence  taken  at  30  frames/second. 
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Table  1 .  Recognition  perfomance 


i 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Our  approach 

17 

20 

21 

21 

21 

22 

22 

23 

24 

24 

Mean-shape  approach 

13 

14 

16 

19 

19 

20 

20 

20 

20 

22 

Landmark-based  approach 

10 

11 

14 

17 

19 

19 

20 

20 

21 

22 

NOTE:  For  each  test  cycle,  we  see  whether  the  correct  class  is  in  the  first  i  classes,  ranked 
according  to  the  metric  in  (3).  This  table  shows  the  number  of  test  cycles  out  of  26  for 
which  that  is  the  case,  plotted  against  i. 


the  mean-shape  approach.  Some  works  in  the  literature  sug¬ 
gest  that  gait  recognition  can  be  achieved  using  merely  a  mean 
shape  of  the  cycle,  rather  than  the  full  cycle  (as  in  Liu  et  al. 
2004b  and  Wang  et  al.  2002).  For  each  person  in  the  training 
set,  these  methods  compute  a  single  mean  shape.  Then,  for  a 
test  sequence,  a  single  mean  shape  is  computed,  and  the  best 
match  in  the  training  set  is  sought.  Surprisingly,  decent  per¬ 
formance  has  been  reported  with  this  simplified  method.  With 
a  slight  variation,  we  computed  two  means  for  each  gait  se¬ 
quence  rather  than  a  single  mean.  We  computed  a  mean  for 
each  person  when  the  arms  and  legs  were  together  and  a  second 
mean  of  the  silhouettes  when  the  arms  and  legs  were  extended. 
Figure  13  shows  a  mean  sequence  of  length  6  as  used  in  the 
previous  method,  along  with  the  two  means  for  the  same  per¬ 
son  that  we  computed  for  the  present  method.  Table  1  summa¬ 
rizes  results  from  this  approach  under  the  mean  shape  method. 
Finally,  we  also  computed  recognition  performance  using  the 
landmark-based  shape  analysis  of  boundary  curves.  Although 
the  general  approach  here  is  same  as  our  method,  the  choices 
of  shape  space  S,  geodesic  lengths  d{-,  ■),  Karcher  means,  and 
so  on  are  different  (Dryden  and  Mardia  1998).  In  this  compu¬ 
tation,  we  represent  each  silhouette  by  100  uniformly  spaced 
pseudolandmarks.  Recognition  results  based  on  this  method  are 
also  reported  in  Table  1 . 

6.  SUMMARY 

In  this  work  we  have  presented  a  novel  framework  for  gait 
recognition,  considering  gait  as  a  cyclostationary  process  on  a 
shape  space  of  simple  closed  curves.  Geometric  tools  allowed 
us  to  perform  interpolation  between  shapes,  registration  of  gait 
cycles,  averaging  of  gait  cycles,  and  comparisons  of  gait  cy¬ 
cles  for  human  recognition.  By  comparing  mean  cycles  rather 
than  the  cycles  themselves,  we  suppressed  intraclass  variability 
and  improved  the  classification  performance.  We  have  demon¬ 
strated  the  classification  technique  on  a  set  of  26  individuals. 


F: jasaap06263rl . tex;  (Aust)  p.  9 

9 

APPENDIX  A:  DIFFERENTIAL 
GEOMETRY:  BACKGROUND 

In  this  section  we  provide  a  short  introduction  to  ideas  from  differ¬ 
ential  geometry  that  are  relevant  to  our  approach.  (For  more  detailed 
discussion,  see,  e.g.,  Spivak  1979;  do  Carmo  1976;  Lang  1999.) 

A  manifold  M  is  a  set  that,  among  other  properties,  has  the  impor¬ 
tant  property  of  being  locally  Euclidean;  that  is,  for  any  point  p  &  M, 
there  exists  a  one-to-one  mapping  between  a  neighborhood  of  p  and  an 
open  subset  of  a  Euclidean  space.  For  any  point  p  €  M,  consider  the 
collection  of  all  parameterized  curves,  c :  (— £,  e)  ^  M,  with  c(0)  =  p. 
Then  TpM  is  the  collection  of  all  vectors  c^(0),  and  it  is  a  vector  space 
even  when  M  is  not  a  vector  space.  If  we  can  define  an  inner  product 
on  TpM  that  varies  smoothly  with  p,  then  M  is  called  a  Riemannian 
manifold,  with  that  inner  product  as  its  Riemannian  metric.  A  group  G 
is  said  to  act  on  M  if  there  exists  a  mapping  G  x  M  M,  denoted  by 
go  M,  that  satisfies  the  following  properties: 

•  We  have  e  o  p  =  p  for  all  p  &  M,  where  e  is  identity  element 
of  G. 

•  If  gi,  g2  e  G,  then  gi  o  (g2  o  p)  =  (gj  ■  g2)  o  p  for  all  p&M. 
Here  gi  •  g2  denotes  the  group  operation  in  G. 

A  distance  measure  on  a  Riemannian  manifold  M  is  indispensable 
for  comparing  elements  of  M.  Because  distances  on  M  are  realized  as 
lengths  of  geodesics,  we  start  by  introducing  geodesics.  A  geodesic  on 
a  Riemannian  manifold  is  a  path  that  is  locally  length- minimizing.  For 
two  points  p,  ^  e  M,  we  let  c :  [0, 1]  ^  M  be  a  differentiable  curve  on 
M  connecting  p  and  q,  that  is,  c(0)  =  p  and  c(l)  =  q.  If  we  let  c'{t) 
be  the  derivative  of  c,  then  the  length  of  the  curve  c  from  p  to  q  is 

lI^{c)  =  j\\c\t)\\gdt, 

where  ||  •  ||g  is  the  norm  induced  by  the  Riemannian  metric  on  M.  The 
distance  between  any  two  points  p  and  q  is  then 

{p,  q)  =  mm{£o(c)  :  c(0)  =  p,  c(l)  =  q}. 

The  curves  c  for  which  this  minimum  is  achieved  are  geodesics  of  the 
manifold  M.  Generally,  geodesics  on  a  M  are  constant-speed  curves  on 
M,  but  this  class  of  curves  may  contain  curves  that  are  not  minimizing. 
For  example,  on  the  manifold  =  {x  e  :  ||x||  =  1},  arcs  of  great 
circles  are  geodesics,  but  a  great-circle  arc  is  minimizing  only  if  its 
length  is  n  or  less.  We  often  denote  4^(/?,  u,  f)  to  be  the  point  obtained 
by  starting  at  p  and  traveling  for  a  time  t  along  a  geodesic  path  in  the 
direction  oi  v  &  TpM. 

On  simpler  manifolds,  such  as  a  sphere,  analytical  expressions  for 
geodesics  are  readily  available.  However,  on  other  manifolds,  includ¬ 
ing  our  shape  space  S,  analytical  derivation  of  a  geodesic  is  intractable, 
and  we  must  use  a  computational  technique.  One  possibility  is  a  shoot¬ 
ing  method,  a  numerical  technique  for  computing  geodesics  on  com¬ 
plicated  nonlinear  manifolds.  Given  two  points  pi,  P2  ^  M,  the  goal 
is  to  compute  a  geodesic  path  between  pi  and  p2.  This  objective  is 


Test  Cycle  Training  Cycle  for  Best  Match 


Test  Cycle  Training  Cycle  for  Same  Person  Training  Cycle  Selected  (Incorrect) 


Figure  12.  Examples  of  a  correct  match  and  an  incorrect  match. 
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Mean  Cycle  of  Length  Six 


The  Pair  of  Mean  Shapes 


Figure  13.  An  example  of  the  training  sequence  used  in  mean-shape-based  classification.  On  the  left  is  a  mean  cycle  of  six  shapes  used  in 
method  1 .  On  the  right  are  the  two  mean  shapes  used  in  the  mean-shape  approach. 


achieved  by  calculating  a  vector  v  e  T'pi  M  such  that  v,  1)  =  P2- 

In  the  shooting  method,  we  begin  with  an  initial  approximation  of  v, 
then  use  a  gradient  method  to  minimize  the  miss  distance  between  suc¬ 
cessive  estimates  of  'Pipi,  v,  1)  and  the  target  point  P2.  The  shooting 
method  for  this  minimization  provides  a  sequence  of  approximations 
to  the  path  that  we  seek.  Figure  1  gives  an  example  of  implementation 
of  the  shooting  method  on  an  ellipsoid,  showing  successive  approxi¬ 
mations  to  the  path  between  two  points. 

A  method  by  which  we  can  go  back  and  forth  between  M  and  Tp  M 
is  the  exponential  map.  This  is  formally  defined  as  follows. 

Definition  2  (Exponential  map).  Let  M  be  a  manifold  and  let  p  e  M 
and  V  €  TpM.  Let  c:  [0, 1]  ^  M  be  the  unique  parameterized  geo¬ 
desic  with  c(0)  =  p  and  c'(0)  =  v.  The  exponential  map  of  u  at  p  is 
expp(u)  =  c(l). 

Definition  3  (Inverse  exponential  map).  Let  M  be  a  manifold  and 
let  p,q  €  M.  We  let  c:  [0,  1]  ^  M  be  a  minimizing  geodesic  with 
c(0)  =  p  and  c(l)  =  ^.  The  inverse  exponential  map  of  q  zl  TpM  is 
any  vector  v  such  that  expp(i;)  =  q. 

We  also  need  the  ability  to  compute  a  sample  mean  shape  on  our 
shape  space.  On  a  manifold,  the  sample  mean  shape  is  taken  to  be  the 
sample  Karcher  mean  defined  as  follows. 

Definition  4  (Karcher  mean).  Suppose  that  Yi, . .  .,Ym  is  a  set  of 
points  on  M.  Then  the  sample  Karcher  mean  set  is  the  set  of  minimiz- 
ers  of  the  function 

^  m 

F(q)  =  -y^d(q,Yk)^,  q€M. 
i=\ 

We  say  that  the  sample  Karcher  mean  exists  if  the  minimizer  is  unique. 

The  Karcher  mean  is  sometimes  called  the  intrinsic  mean,  because 
it  uses  an  intrinsic  metric  on  a  manifold  as  opposed  to  an  extrinsic 
mean  computed  by  embedding  a  manifold  in  an  ambient  Euclidean 
space.  There  exist  iterative  techniques  for  finding  the  Karcher  mean 
that  use  the  gradient  of  F  to  update  the  estimate  until  convergence.  A 
commonly  used  method  updates  the  mean  /x  as 

y  m 

ji  exp^(u),  where  v  =  -  ^exp^HL). 

/  =  1 

APPENDIX  B:  DYNAMIC  PROGRAMMING 

Dynamic  programming  is  an  important  tool  in  the  registration  of 
curves  and  sequences,  and  it  can  help  solve  the  following  general  prob¬ 
lem.  For  two  given  functions  / :  [0,  a]  i-^  K  and  h :  [0,  b]  K,  find  a 
function  g :  [0,  a]  [0,  b]  that  minimizes  the  functional, 

{f{x)-h{g(x))f  dx. 

In  addition,  g  is  assumed  to  be  differentiable  and  to  have  a  positive 
derivative  at  each  point.  An  overview  of  the  DP  algorithm  is  as  follows. 


We  divide  the  region  [0,  a]  x  [0,  b\  into  a  grid  with  N  x  N  nodes.  Ex¬ 
plicitly,  the  node  (/,  j)  in  the  grid  represents  the  point 
We  obtain  a  piecewise-linear  approximation  to  g  on  [0,  a]  and  con- 
stmct  a  g  that  passes  through  nodes  on  the  grid. 

Eor  any  node  in  the  grid,  define  the  cost  H{i,  j)  at  node  (i,  j)  to  be 
the  minimum  cost  of  reaching  (i,  j)  on  a  strictly  increasing  path  from 
(0, 0),  and  set  H(0,  0)  =  0.  Now  our  problem  is  to  find  a  path  that 
achieves  H(N  —  1,  A  —  1).  To  minimize  the  cost,  H(i,  j),  of  reaching 
(i,  j),  we  first  denote  A(/ y)  to  be  the  set  of  nodes  that  are  permissible 
predecessors  to  the  node  (/,  j).  Elements  of  N(i  j)  are  selected  to  en¬ 
sure  that  the  slope  of  a  candidate  path  g  is  positive  at  each  point  and 
are  close  to  (;,  j). 

Eor  (k,l)  €  A(,' j),  the  cost  of  traveling  to  (;, /)  through  a  path 
that  contains  {k,  1)  is  the  minimum  cost  of  traveling  to  (k,  1),  plus 
the  cost  of  traveling  from  (k,/)  to  {i,i)  along  a  linear  path.  Let¬ 
ting  E((i,j),  {k,l))  denote  the  cost  of  traveling  to  (/,  7)  through  a 
path  that  contains  (k,l),  this  says  that  E((i,  j),{k,l))  =  H(k,l)  + 
L{{k,  1),  (i,  j)),  where  L{{k,  1),  (i,  j))  is  the  cost  of  traveling  from 
(k,  1)  to  ((,  j)  in  a  straight  line  connecting  those  two  points.  Now  the 
cost  of  reaching  node  (/,  j)  is  F[{k,  /),  where 

(k,  /)  =  argmin  E{{k,  /),  (;,  j)). 

We  find  g  by  computing  //(n  —  1, «  —  1)  in  an  iterative  manner.  In 
doing  so,  we  obtain  nodes  (mj,  ui), . . . ,  V]fi)  =  {N  —  I,  N  —  1), 
representing  the  critical  path. 

[Received  May  2006.  Revised  January  2007.] 
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